//
// nnMeteo.cc
//
// programma di esempio per l'uso e la verifica dei risultati di una
// rete neurale addestrata a fare previsioni metereologiche
//
// Istruzioni per l'uso:
//
// 1) eseguire root:
// shell> root <return>
//
// 2) compilare ed eseguire il programma da dentro Root
// root [0] .x nnMeteo.cc+
//
// (per eseguire in una volta sola la conversione dei dati (da ascii a
// root-tree), la visualizzazione delle variabili, e l'addestramento e
// utilizzo della rete.)
//
// 3) oppure, sempre da dentro Root:
//
// root [1] .L nnMeteo.cc+              (per caricare in memoria la macro; farlo tutte e sole le volte che la si e` modificata)
// root [2] converti("pippo",K)         (per convertire pippo.txt in pippo.root, su K giorni)
// root [3] visualizza("pippo")         (per visualizzare le variabili ora contenute in pippo.root)
// root [4] reteneurale("pippo",K,N)    (per addestrare la rete sui dati di pippo.root con N cicli di addestramento, con l'algoritmo K, e salvarla nel file pesi.txt)
// root [5] usarete("pippo",K,N)        (per utilizzare la rete gia' addestrata di cui avete salvato i pesi in pesi.txt, sul sample pippo.root)
//

#if !defined(__CINT__) || defined(__MAKECINT__)

#include <assert.h>
#include <fstream>
#include <iostream>
#include <iomanip>
#include <string>
#include <vector>
#include <numeric>

#include <string.h>

#include "TROOT.h"
#include "TSystem.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TLegend.h"
#include "TH1F.h"
#include "TMath.h"
#include "TFile.h"
#include "TTree.h"
#include "TEventList.h"
#include "TMultiLayerPerceptron.h"
#include "TMLPAnalyzer.h"

#endif

using namespace std;

//
// meteoData: classe per oggetto contenente i dati metereologici
//
// le macro ClassDef e ClassImp servono perche' Root possa
// usare le variabili e i metodi della classe per srivere/leggere
// su disco e per le formule [ es. ntupla->Draw() ]
//
// i dati si possono accedere sia con il vettore data[]
// sia con i metodi
// - pressione()
// - temperatura()
// - umidita()
// - precipitazione()
//
// le parentesi () sono necessarie ma a parte questo questi metodi si
// usano come se fossero variabili, cioe' si possono usare per fare
// assegnamenti del tipo "md.pressione() = 1000;"
//
class meteoData {
 private:
 public:
  enum { meteoVarNum = 4 };
  Double_t data[meteoVarNum];

  //--- access for reading and writing
  Double_t& pressione() {return data[0];}
  Double_t& temperatura() {return data[1];}
  Double_t& umidita() {return data[2];}
  Double_t& precipitazione() {return data[3];}

  //--- access for reading when object is declared constant
  Double_t const& pressione() const {return data[0];}
  Double_t const& temperatura() const {return data[1];}
  Double_t const& umidita() const {return data[2];}
  Double_t const& precipitazione() const {return data[3];}

  meteoData() { memset( &data[0], sizeof(data), 0); }
  virtual ~meteoData() {}
  ClassDef(meteoData,1)  // per segnalare a Root il nome delle variabili
};

ClassImp(meteoData);

//
// salva i dati raccolti nel vettore in un file di Root
//
void salvaNupla(TString name, vector<meteoData>& mdl,
		Int_t nSkip, Int_t nSkipMax )
{
  //--- il nome del file di Root dipende da quale giorno si parte
  TString output = name;
  //--- nome dell'oggetto nella n-upla
  TString mdName = "md";

  if (nSkip != 0) {
    output += "-";
    output += nSkip;
    mdName += nSkip;
  }
  output += ".root";
  
  // ora dobbiamo salvare tutto in un "tree" di root
  // apertura file di output
  TFile rfile(output,"RECREATE");

  TTree tree("T", "il nostro tree");
  meteoData* mdp = new meteoData();
  
  tree.Branch(mdName, "meteoData", &mdp, 64000, 1);
  
  //
  // produciamo un file contenente i dati corrispondenti ai giorni
  // da <nSkip> a <ndays - (nSKipMax-nSkip) -1> (inclusi)
  //
  Int_t ndays = mdl.size();
  Int_t iDay = 0;
  Int_t iValidDay = 0;
  
  vector<meteoData>::iterator iter;
  for( iter = mdl.begin(); iter != mdl.end(); iter++, iDay++ ) {
    //
    // salviamo nell'oggetto puntato da mdp
    // l'-i-esimo giorno di dati metereologici
    //
    *mdp = *iter;
    
    //
    // eliminiamo i dati non validi
    // questo significa che in alcuni casi il "giorno precedente"
    // corrispondera' in realta' a diversi giorni prima
    // non e' facile evitare questo senza significative
    // modifiche al programma
    //
    if (mdp->pressione() == -9999 ||
	mdp->temperatura() == -9999 ||
	mdp->umidita() == -9999 ||
	mdp->precipitazione() == -9999) {
      continue;
    }
    
    //--- qui si possono regolarizzare i dati

    // mdp->precipitazione() = mdp->precipitazione() != 0 ? 1 : 0;
    
    //--- contatore per i soli giorni con dati validi
    if (iValidDay++ < nSkipMax -nSkip) continue;
    if (iDay >= ndays - nSkip) continue;
    //--- mettiamo i valori nella n-upla
    tree.Fill();
    // cout << "day " << iDay << " pressione " << mdp->pressione() << endl;
  }
  
  // operazioni conclusive: chiusura file di input, scrittura degli output
  tree.Write();
  rfile.Write();
  rfile.Close();
  delete mdp;
  
  cout << "Scritto il Root file " << output << endl;
}

//
// converte il file di testo name.txt nel file di Root name.root
// attenzione: se invocato piu' volte il file di Root risultante
// corrispondera' ai dati del solo ultimo file di dati
//
void converti(TString name, Int_t ndays)
{
  // definizione dei nomi dei file di input e output
  TString input = name + ".txt";

  // apertura file dati per leggere
  ifstream in;
  in.open(input);

  // lettura del contenuto del file di input
  vector<meteoData> mdl(ndays);
  meteoData md;
  
  Int_t nlines = 0;
  Double_t dummy;

  for(;;) {
    for (Int_t i=0; i<ndays; i++) {
      //
      // ogni linea corrisponde a 365 dati di una variabile
      // --> si aggiorna il campo corrispondente nell' i-esimo elemento
      // del vettore, la corrispondenza e':
      // linea 1 -> pressione
      // linea 2 -> temperatura
      // linea 3 -> umidita'
      // linea 4 -> precipitazione
      //
      if (nlines < meteoData::meteoVarNum) {
	in >> mdl[i].data[nlines];
      } else {
	// --- dati non utilizzati nelle linee successive a meteoVarNum
	in >> dummy;
      }
    }
    if (!in.good()) break;
    nlines++;
  }
  in.close();

  cout << "Letti dati corrispondenti a " << nlines << " giorni" << endl;

#if 0
  vector<meteoData>::iterator iter;
  for( iter = mdl.begin(); iter != mdl.end(); iter++) {
    cout
      << "p, t, h, r = "
      << iter->pressione() << " "
      << iter->temperatura() << " "
      << iter->umidita() << " "
      << iter->precipitazione() << " "
      << endl;
  }
#endif

  //
  // salva i dati in diversi file di Root contenenti gli stessi valori
  // ma sfasati da 0 a 4 giorni
  //
  salvaNupla(name, mdl, 0, 4);
  salvaNupla(name, mdl, 1, 4);
  salvaNupla(name, mdl, 2, 4);
  salvaNupla(name, mdl, 3, 4);
  salvaNupla(name, mdl, 4, 4);
  
  return;
}

void visualizza(TString name)
{
  TString rfname = name + ".root";

  // apertura del file e del suo tree
  TFile rfile(rfname);  
  // metto "T" perche' e' il nome del tree contenuto nel file
  TTree *tree = (TTree*)gDirectory->Get("T");

  TCanvas* c1 = new TCanvas("c1","Variabili",40,40,800,800);
  c1->SetFillColor(18);
  TPad* pad1 = new TPad("pad1","Pressione",0.05,0.50,0.45,0.95,21);
  TPad* pad2 = new TPad("pad2","Temperatura",0.05,0.05,0.45,0.45,21);
  TPad* pad3 = new TPad("pad3","Umidita`",0.50,0.50,0.95,0.95,21);
  TPad* pad4 = new TPad("pad4","Precipitazione",0.50,0.05,0.95,0.45,21);
  pad1->Draw();
  pad2->Draw();
  pad3->Draw();
  pad4->Draw();
  pad1->cd();
  tree->Draw("pressione()");
  pad2->cd();
  tree->Draw("temperatura()");
  pad3->cd();
  tree->Draw("umidita()");
  pad4->cd();
  tree->Draw("precipitazione()");
}

TString nomeAlgo(Int_t nAlgo)
{
  TString algo;
  if (nAlgo == 1) algo="Stochastic";
  if (nAlgo == 2) algo="Batch";
  if (nAlgo == 3) algo="SteepestDescent";
  if (nAlgo == 4) algo="RibierePolak";
  if (nAlgo == 5) algo="FletcherReeves";
  if (nAlgo == 6) algo="BFGS";

  return algo;
}

TString nomeBase(TString name,
		 Int_t nAlgo, Int_t neuro1, Int_t neuro2, Int_t ntrain)
{
  // salva i pesi di questa rete per poterli usare dopo:
  TString base = name + "_";
  base += nomeAlgo(nAlgo);
  base += "_";
  base += neuro1;
  base += "_";
  base += neuro2;
  base += "_";
  base += ntrain;

  return base;
}

TString nomeFilepesi(TString name,
		     Int_t nAlgo, Int_t neuro1, Int_t neuro2,Int_t ntrain)
{
  return
    nomeBase(name, nAlgo, neuro1, neuro2, ntrain)
    + "_pesi.txt";
}

TString nomeNnConfig(TString name,
		     Int_t nAlgo, Int_t neuro1, Int_t neuro2,Int_t ntrain)
{
  return
    nomeBase(name, nAlgo, neuro1, neuro2, ntrain)
    + "_config.txt";
}

//
// costruisci la stringa che definisce la rete neurale
//
TString nnStruct(Int_t neuro1, Int_t neuro2)
{
  // istruzioni qui:
  // http://root.cern.ch/root/html310/TMultiLayerPerceptron.html

  TString struttura("");

//   struttura += "md1.pressione()-md2.pressione(),";
//   struttura += "md1.temperatura()-md2.temperatura(),";
//   struttura += "md1.umidita()-md2.umidita(),";

  struttura += "md1.pressione(),";
  struttura += "md1.temperatura(),";
  struttura += "md1.umidita(),";
  struttura += "md1.precipitazione()";

  struttura += ",";

  struttura += "md2.pressione(),";
  struttura += "md2.temperatura(),";
  struttura += "md2.umidita(),";
  struttura += "md2.precipitazione()";

//   struttura += ",";
//
//   struttura += "md3.pressione(),";
//   struttura += "md3.temperatura(),";
//   struttura += "md3.umidita(),";
//   struttura += "md3.precipitazione()";
//
//   struttura += ",";
//
//   struttura += "md4.pressione(),";
//   struttura += "md4.temperatura(),";
//   struttura += "md4.umidita(),";
//   struttura += "md4.precipitazione()";

  struttura += ":";
  struttura += neuro1;
  if (neuro2>0) { // due strati nascosti
    struttura+=":";
    struttura+=neuro2;
  }
  struttura+=":";
  //
  // output: la @ significa che la rete deve normalizzarlo
  // gli input sono sempre automaticamente normalizzati
  //
  struttura+="@md.precipitazione()";
  
  return struttura;
}

void reteneurale(TString name,
		 Int_t nAlgo, Int_t neuro1, Int_t neuro2, Int_t ntrain)
{
  if (!gROOT->GetClass("TMultiLayerPerceptron")) {
    gSystem->Load("libMLP");
    cout
      << endl
      << "E' stata caricata la libreria ROOT per le reti neurali."
      << endl
      << endl;
  }
  
  //
  // apertura del file e del suo tree
  // questo sara' il file usato per il training della rete
  //
  TFile rfile(name + ".root");
  // metto "T" perche' e' il nome del tree contenuto nel file
  TTree *tree = (TTree*)rfile.Get("T");

  //
  // apriamo i file di Root sfasati di 1-4 giorni
  // e usiamoli come file "amici" del precedente,
  // per poter usare nella rete neurale i dati dei giorni
  // precedenti
  //
  tree->AddFriend("T", name + "-1.root");
  tree->AddFriend("T", name + "-2.root");
  tree->AddFriend("T", name + "-3.root");
  tree->AddFriend("T", name + "-4.root");

  // creiamo un oggetto della classe "rete neurale" istruzioni qui:
  // http://root.cern.ch/root/html310/TMultiLayerPerceptron.html
  TString struttura =  nnStruct(neuro1, neuro2);
  TMultiLayerPerceptron *mlp = new TMultiLayerPerceptron(struttura, tree);

  TEventList train;
  TEventList test;
  for (int i = 0; i < tree->GetEntries(); i++) {
    if (i%16 < 8)
      train.Enter(i);
    else
      test.Enter(i);
  }
  train.Print();
  test.Print();

  mlp->SetTrainingDataSet(&train);
  mlp->SetTestDataSet(&test);

  //--- algoritmo di addestramento:
  if (nAlgo == 1) mlp->SetLearningMethod(TMultiLayerPerceptron::kStochastic);
  if (nAlgo == 2) mlp->SetLearningMethod(TMultiLayerPerceptron::kBatch);
  if (nAlgo == 3) mlp->SetLearningMethod(TMultiLayerPerceptron::kSteepestDescent);
  if (nAlgo == 4) mlp->SetLearningMethod(TMultiLayerPerceptron::kRibierePolak);
  if (nAlgo == 5) mlp->SetLearningMethod(TMultiLayerPerceptron::kFletcherReeves);
  if (nAlgo == 6) mlp->SetLearningMethod(TMultiLayerPerceptron::kBFGS);
  if (nAlgo<1 || nAlgo>6) cout << "Algoritmo non valido." << endl;
  cout << "Algoritmo scelto: " << nomeAlgo(nAlgo) << endl;

  // addestramento:
  mlp->Randomize(); // randomizza i pesi tra neuroni prima dell'addestramento
  mlp->Train(ntrain, "text,graph,update=10");

  // Use TMLPAnalyzer to see what it looks for
  TCanvas* mlpa_canvas = new TCanvas("mlpa_canvas", "Network analysis");
  mlpa_canvas->SetWindowPosition(160,80);
  TMLPAnalyzer ana(mlp);
  // Initialisation
  ana.GatherInformations();
  // output to the console
  ana.CheckNetwork();
  mlp->Draw();

  // mlp->DrawResult(0, "test"); // perche' non funziona???

  TString nomepesi = nomeFilepesi(name, nAlgo, neuro1, neuro2, ntrain);
  mlp->DumpWeights(nomepesi);
  cout << "File dei pesi salvato col nome " << nomepesi << endl;
  // per scoprire come fissare l'algoritmo di addestramento, e in
  // generale come usare la rete, leggere attentamente
  // http://root.cern.ch/root/html310/TMultiLayerPerceptron.html
}

void usarete(TString name,
	     Int_t nAlgo, Int_t neuro1, Int_t neuro2, Int_t ntrain)
{
  if (!gROOT->GetClass("TMultiLayerPerceptron")) {
    gSystem->Load("libMLP");
    cout
      << endl
      << "E' stata caricata la libreria ROOT per le reti neurali."
      << endl
      << endl;
  }

  //
  // apertura del file e del suo tree
  // questo sara' il file usato per il training della rete
  //
  TFile rfile(name + ".root");
  // metto "T" perche' e' il nome del tree contenuto nel file
  TTree *tree = (TTree*)rfile.Get("T");

  //
  // apriamo i file di Root sfasati di 1-4 giorni
  // e usiamoli come file "amici" del precedente,
  // per poter usare nella rete neurale i dati dei giorni
  // precedenti
  //
  tree->AddFriend("T", name + "-1.root");
  tree->AddFriend("T", name + "-2.root");
  tree->AddFriend("T", name + "-3.root");
  tree->AddFriend("T", name + "-4.root");

  meteoData* mdp = 0;
  meteoData* mdp1 = 0;
  meteoData* mdp2 = 0;
  meteoData* mdp3 = 0;
  meteoData* mdp4 = 0;
  tree->SetBranchAddress("md", &mdp);
  tree->SetBranchAddress("md1", &mdp1);
  tree->SetBranchAddress("md2", &mdp2);
  tree->SetBranchAddress("md3", &mdp3);
  tree->SetBranchAddress("md4", &mdp4);

  // definiamo la struttura della rete neurale (dev'essere identica a
  // quella usata in reteneurale())
  TString struttura =  nnStruct(neuro1, neuro2);
  TMultiLayerPerceptron *mlp =
    new TMultiLayerPerceptron(struttura, tree);

  mlp->LoadWeights(nomeFilepesi(name, nAlgo, neuro1, neuro2, ntrain));
  //  mlp->Draw(); // questo in realta' non ci serve

  // loop sul numero di entries nel sample
  Double_t inputs[4*4];
  Double_t output;
  Double_t target;
  TH1F *h1;
  TH1F *h2;
  h1 = new TH1F("h1","output",20,-1,10);
  h2 = new TH1F("h2","target",20,-1,10);
  Int_t pioveNN=0;
  Int_t pioveR=0;
  Int_t falsiPos=0;
  Int_t falsiNeg=0;

  // soglia sotto la quale non consideriamo piovosa la giornata
  Double_t soglia=1;
  cout << "Numero giorni totali: " << tree->GetEntries() << endl;
  for (int i = 0; i < tree->GetEntries(); i++) {
    tree->GetEntry(i);

    Int_t ii(0);

//     inputs[ii++] = mdp1->pressione()-mdp2->pressione();
//     inputs[ii++] = mdp1->temperatura()-mdp2->temperatura();
//     inputs[ii++] = mdp1->umidita()-mdp2->umidita();

    inputs[ii++] = mdp1->pressione();
    inputs[ii++] = mdp1->temperatura();
    inputs[ii++] = mdp1->umidita();
    inputs[ii++] = mdp1->precipitazione();
    
    inputs[ii++] = mdp2->pressione();
    inputs[ii++] = mdp2->temperatura();
    inputs[ii++] = mdp2->umidita();
    inputs[ii++] = mdp2->precipitazione();

//    inputs[ii++] = mdp3->pressione();
//    inputs[ii++] = mdp3->temperatura();
//    inputs[ii++] = mdp3->umidita();
//    inputs[ii++] = mdp3->precipitazione();
//
//    inputs[ii++] = mdp4->pressione();
//    inputs[ii++] = mdp4->temperatura();
//    inputs[ii++] = mdp4->umidita();
//    inputs[ii++] = mdp4->precipitazione();

    output=mlp->Evaluate(0,inputs);
    target = mdp->precipitazione();

    h1->Fill(target);
    h2->Fill(output);

    if (output > soglia) pioveNN++;
    if (target > soglia) pioveR++;
    if (output > soglia && target < soglia) falsiPos++;
    if (output < soglia && target > soglia) falsiNeg++;
  }

  // visualizza confronto tra output e target
  gStyle->SetOptStat(0000);
  TCanvas* c = new TCanvas("c","Output / Target comparison"); 
  c->SetWindowPosition(260,100);
  c->cd();
  h1->SetLineWidth(2);
  h2->SetLineWidth(2);
  h1->SetLineColor(kBlack);
  h2->SetLineColor(kRed);
  h2->DrawCopy();
  h1->DrawCopy("same");
  //legenda
  Double_t ymax = TMath::Max(h1->GetMaximum(),h2->GetMaximum());
  TLegend *legend =
    new TLegend(7.5,0.75*ymax,9.5,0.95*ymax,"Precipitazioni","br");
  legend->AddEntry(h1,"previsione","l");
  legend->AddEntry(h2,"realta`","l");
  legend->Draw();

  cout
    << "Ha piovuto in " << pioveR
    << " giorni contro " << pioveNN
    << " previsti dalla rete."
    << endl
    << "Falsi positivi: " << falsiPos
    << ", falsi negativi: " << falsiNeg
    << endl;
}

void nnMeteo()
{
  gROOT->Reset(); // fa il reset di tutte le variabili

  // 
  // converte il file testo in un file di ROOT; la variabile numerica
  // e' il numero di giorni descritti dal file se avete piu' di un
  // file di input, ripetete il rigo sopra piu' volte mettendo ogni
  // volta il nome del file (senza ".txt")
  converti("dati", 365);
  
  // grafico delle distribuzioni di tutte le variabili
  visualizza("dati");

  //
  // argomenti di reteneurale()
  // - nome file Root
  // - algoritmo
  //     1) Stochastic
  //     2) Batch
  //     3) SteepestDescent
  //     4) RibierePolak
  //     5) FletcherReeves
  //     6) BFGS
  // - numero neuroni primo strato
  // - numero neuroni secondo strato (puo' essere zero)
  // - numero di iterazioni
  //

  // addestra la rete con il primo algoritmo, 4 neuroni nel primo
  // strato, 2 nel secondo strato, e dopo l'addestramento salva i pesi
  // in un file di testo
  reteneurale("dati",6,2,0,50);

  // ricrea la rete con i pesi dati dall'addestramento, e la mette a
  // confronto con la realta'
  usarete("dati",6,2,0,50);
}
